## this script does the following:
	## 1) combines the model estimates and calculates the standard errors that appear in Table 2 of the article ##
	## 2) uses the posterior draws from the full models to produce the plots in Figure 1 ##
	## 3) aggregates the simulated vote losses generated in scripts "sim" 1-10 to produce the Figure 2 plot##
## read in data ##
## get that data in order ##

	rm(list = ls())
	library(foreign)
	setwd('replication')
	all <- read.table('results.csv', sep = ',', header = TRUE, skip = 1)
	all <- all[, -1]
	all <- all[-15, ]

## now calculate those paramters and standard errors ##
## simple model ##
	
	data <- all[, 1:10]
	m <- 10
	coef <- rowMeans(data[seq(1, 9, 2), ])
	sd <- rowMeans(abs(data[seq(15, 23, 2), ]))
	
## record variance in parameter estimates across imputations ##
	s2q1 <- sum((coef[1] - data[1, ])^2) / (m - 1)
	s2q2 <- sum((coef[2] - data[3, ])^2) / (m - 1)
	s2q3 <- sum((coef[3] - data[5, ])^2) / (m - 1)
	s2q4 <- sum((coef[4] - data[7, ])^2) / (m - 1)
	s2q5 <- sum((coef[5] - data[9, ])^2) / (m - 1)

	s2q1sd <- sum((sd[1] - data[15, ])^2) / (m - 1)
	s2q2sd <- sum((sd[2] - data[17, ])^2) / (m - 1)
	s2q3sd <- sum((sd[3] - data[19, ])^2) / (m - 1)
	s2q4sd <- sum((sd[4] - data[21, ])^2) / (m - 1)
	s2q5sd <- sum((sd[5] - data[23, ])^2) / (m - 1)
	
## calculate standard errors leveraging individual standard errors ##
## and cross-imputation parameter variance ##
	se1 <- sqrt(mean(data[2, ]^2) + s2q1*(1+(1/m)))
	se2 <- sqrt(mean(data[4, ]^2) + s2q2*(1+(1/m)))
	se3 <- sqrt(mean(data[6, ]^2) + s2q3*(1+(1/m)))
	se4 <- sqrt(mean(data[8, ]^2) + s2q4*(1+(1/m)))
	se5 <- sqrt(mean(data[10, ]^2) + s2q5*(1+(1/m)))

	se1sd <- sqrt(mean(data[16, ]^2) + s2q1sd*(1+(1/m)))
	se2sd <- sqrt(mean(data[18, ]^2) + s2q2sd*(1+(1/m)))
	se3sd <- sqrt(mean(data[20, ]^2) + s2q3sd*(1+(1/m)))
	se4sd <- sqrt(mean(data[22, ]^2) + s2q4sd*(1+(1/m)))
	se5sd <- sqrt(mean(data[24, ]^2) + s2q5sd*(1+(1/m)))

	se <- c(se1, se2, se3, se4, se5)
	sdSe <- c(se1sd, se2sd, se3sd, se4sd, se5sd)

## assign data ##
	names <- c('gov', 'govCompSub', 'lr', 'lrGov', 'economyGov')
	
	simple <- cbind(names, round(coef, 3), round(se, 3), abs(round(sd, 3)), round(sdSe, 3))
	simple

## now the full model ##
	
	data <- all[, 11:20]
	m <- 10
	coef <- rowMeans(data[seq(1, 13, 2), ])
	sd <- rowMeans(abs(data[seq(15, 27, 2), ]))
	
## record variance in parameter estimates across imputations ##
	s2q1 <- sum((coef[1] - data[1, ])^2) / (m - 1)
	s2q2 <- sum((coef[2] - data[3, ])^2) / (m - 1)
	s2q3 <- sum((coef[3] - data[5, ])^2) / (m - 1)
	s2q4 <- sum((coef[4] - data[7, ])^2) / (m - 1)
	s2q5 <- sum((coef[5] - data[9, ])^2) / (m - 1)
	s2q6 <- sum((coef[6] - data[11, ])^2) / (m - 1)
	s2q7 <- sum((coef[7] - data[13, ])^2) / (m - 1)

	s2q1sd <- sum((sd[1] - data[15, ])^2) / (m - 1)
	s2q2sd <- sum((sd[2] - data[17, ])^2) / (m - 1)
	s2q3sd <- sum((sd[3] - data[19, ])^2) / (m - 1)
	s2q4sd <- sum((sd[4] - data[21, ])^2) / (m - 1)
	s2q5sd <- sum((sd[5] - data[23, ])^2) / (m - 1)
	s2q6sd <- sum((sd[6] - data[25, ])^2) / (m - 1)
	s2q7sd <- sum((sd[7] - data[27, ])^2) / (m - 1)
	
## calculate standard errors leveraging individual standard errors ##
## and cross-imputation parameter variance ##
	se1 <- sqrt(mean(data[2, ]^2) + s2q1*(1+(1/m)))
	se2 <- sqrt(mean(data[4, ]^2) + s2q2*(1+(1/m)))
	se3 <- sqrt(mean(data[6, ]^2) + s2q3*(1+(1/m)))
	se4 <- sqrt(mean(data[8, ]^2) + s2q4*(1+(1/m)))
	se5 <- sqrt(mean(data[10, ]^2) + s2q5*(1+(1/m)))
	se6 <- sqrt(mean(data[12, ]^2) + s2q6*(1+(1/m)))
	se7 <- sqrt(mean(data[14, ]^2) + s2q7*(1+(1/m)))

	se1sd <- sqrt(mean(data[16, ]^2) + s2q1sd*(1+(1/m)))
	se2sd <- sqrt(mean(data[18, ]^2) + s2q2sd*(1+(1/m)))
	se3sd <- sqrt(mean(data[20, ]^2) + s2q3sd*(1+(1/m)))
	se4sd <- sqrt(mean(data[22, ]^2) + s2q4sd*(1+(1/m)))
	se5sd <- sqrt(mean(data[24, ]^2) + s2q5sd*(1+(1/m)))
	se6sd <- sqrt(mean(data[26, ]^2) + s2q6sd*(1+(1/m)))
	se7sd <- sqrt(mean(data[28, ]^2) + s2q7sd*(1+(1/m)))

	se <- c(se1, se2, se3, se4, se5, se6, se7)
	sdSe <- c(se1sd, se2sd, se3sd, se4sd, se5sd, se6sd, se7sd)

## assign data ##
	names <- c('gov', 'govCompSub', 'lr', 'lrGov', 'economyGov', 'compLr', 'comEcon')
	
	full <- cbind(names, round(coef, 3), round(se, 3), abs(round(sd, 3)), round(sdSe, 3))
	full

################################################################################################

## create graph of effects of compromise on left-right weights ##
## this is for germany ##
	bbt <- c()
	
	for(i in 1:10){
		if(i == 1){
			bbt <- na.omit(read.dta(paste('fullBeta', i, '.dta', sep = '')))			
			data <- read.dta(paste('mixBetas', i, '.dta', sep = ''))
			data <- data[, -c(38:41)]
			data$iteration <- i
		}else{
			bbt <- rbind(na.omit(read.dta(paste('fullBeta', i, '.dta', sep = ''))), bbt)
			dataHold <- read.dta(paste('mixBetas', i, '.dta', sep = ''))
			dataHold <- dataHold[, -c(38:41)]
			dataHold$iteration <- i
			data <- rbind(data, dataHold)
		}
	}

	bbt <- bbt[, c(1:7)]

	pdf('marginalEffects.pdf', width = 14)
	par(mfrow = c(1, 2), oma = c(0, 0, 3, 0))
	## marginal effects for economy ###

	range <- seq(0, 6, by = 0.025)
	
	supHi <- c()
	supMed <- c()
	supLo <- c()

	compBase <- mean(data$govCompSub[data$gov == 1])
	compNew <- mean(data$govCompSub[data$gov == 1]) + sd(data$govCompSub[data$gov == 1])
	data$economy <- data$economy + abs(min(data$economy))
	econ <- mean(data$economy[data$gov == 1])
	lr <- mean(data$lr)

	for(i in 1:length(range)){

		supGov <- exp(bbt[, 1] + 
			bbt[, 2] * compBase +
			bbt[, 3] * lr +
			bbt[, 4] * lr + 
			bbt[, 5] * range[i] + 
			bbt[, 6] * compBase * lr +
			bbt[, 7] * compBase * range[i])

		supGovNew <- exp(bbt[, 1] + 
			bbt[, 2] * compNew +
			bbt[, 3] * lr +
			bbt[, 4] * lr + 
			bbt[, 5] * range[i] + 
			bbt[, 6] * compNew * lr +
			bbt[, 7] * compNew * range[i])


		supOpp <- exp(bbt[, 3] * lr)

		supProb <- supGov/(supGov + supOpp)
		supProbNew <- supGovNew/(supGovNew + supOpp)
		supEffect <- supProbNew - supProb

		supHi[i] <- quantile(supEffect, 0.975)
		supMed[i] <- quantile(supEffect, 0.5)
		supLo[i] <- quantile(supEffect, 0.025)

	}
		
		plot(1, 1, type = 'n',
			xlim = c(min(range), max(range)),
			ylim = c(-0.15, 0.1),
			ylab = '',
			xlab = 'Observed Range of Economic Performance',
			main = '',
			axes = FALSE)
		polygon(c(0, 6, 6, 0), c(0.1, 0.1, -0.15, -0.15), col = rgb(0, 0, 0, 0.05), border = FALSE)
		for(i in 1:12){
			lines(c(i/2, i/2), c(0.1, -0.15), col = 'white', lwd = 2)
		}
		for(i in 1:10){
			lines(c(0, 10), c(-0.15+0.025*i, -0.15+0.025*i), col = 'white', lwd = 2)
		}
		lines(c(0, 6), c(0, 0), col = rgb(0, 0, 0, 0.5), lwd = 2, lty = 3)
		drop <- data$economy
		x <- seq(0, 6, by = 0.5)
		xCount <- c()
		for(i in 1:(length(x) - 1)){
			xCount[i] <- length(drop[drop >= x[i] & drop < x[i+1]])
		}
		for(i in 1:(length(x) - 1)){
			polygon(c(x[i], x[i+1], x[i+1], x[i]), c(-0.15, -0.15, -0.15+xCount[i]/max(xCount/0.25), -0.15+xCount[i]/max(xCount/0.25)), col = rgb(0, 0, 0, 0.1), border = FALSE)
		}
		axis(1)
		axis(2, line = -1)
		axis(4, at = c(-0.15, 0.1), labels = c('Min', 'Max'), line = -1)
		polygon(c(range, rev(range)), c(supHi, rev(supLo)), col = rgb(0, 0, 0, 0.25), border = FALSE)
		lines(range, supMed, col = rgb(0, 0, 0, 1), lwd = 2)
		mtext('Marginal Effect of Perceived Compromise on Incumbent Selection\nOver Economic Performance', font = 2)
		mtext('Observed Values of Economic Performance', side = 4)
		mtext('Change in Probability of Choosing Incumbent Party', side = 2, line = 1.5)

	range <- seq(0, 5, by = 0.025)
	
	supHi <- c()
	supMed <- c()
	supLo <- c()

	compBase <- mean(data$govCompSub[data$gov == 1])
	compNew <- mean(data$govCompSub[data$gov == 1]) + sd(data$govCompSub[data$gov == 1])
	econ <- mean(data$economy[data$gov == 1])

	for(i in 1:length(range)){

		supGov <- exp(bbt[, 1] + 
			bbt[, 2] * compBase +
			bbt[, 3] * range[i] +
			bbt[, 4] * range[i] + 
			bbt[, 5] * econ + 
			bbt[, 6] * compBase * range[i] +
			bbt[, 7] * compBase * econ)

		supGovNew <- exp(bbt[, 1] + 
			bbt[, 2] * compNew +
			bbt[, 3] * range[i] +
			bbt[, 4] * range[i] + 
			bbt[, 5] * econ + 
			bbt[, 6] * compNew * range[i] +
			bbt[, 7] * compNew * econ)


		supOpp <- exp(bbt[, 3] * range[i])

		supProb <- supGov/(supGov + supOpp)
		supProbNew <- supGovNew/(supGovNew + supOpp)
		supEffect <- supProbNew - supProb

		supHi[i] <- quantile(supEffect, 0.975)
		supMed[i] <- quantile(supEffect, 0.5)
		supLo[i] <- quantile(supEffect, 0.025)

	}

		plot(1, 1, type = 'n',
			xlim = c(min(range), max(range)),
			ylim = c(-0.2, 0.05),
			ylab = '',
			xlab = 'Observed Range of Ideological Distance Between Voters and Parties',
			main = '',
			axes = FALSE)
		polygon(c(0, 5, 5, 0), c(0.05, 0.05, -0.2, -0.2), col = rgb(0, 0, 0, 0.05), border = FALSE)
		for(i in 1:10){
			lines(c(i/2, i/2), c(0.05, -0.2), col = 'white', lwd = 2)
		}
		for(i in 1:12){
			lines(c(0, 10), c(-0.2+0.025*i, -0.2+0.025*i), col = 'white', lwd = 2)
		}
		lines(c(0, 5), c(0, 0), col = rgb(0, 0, 0, 0.5), lwd = 2, lty = 3)
		drop <- data$lr
		x <- seq(0, 5, by = 0.25)
		xCount <- c()
		for(i in 1:(length(x) - 1)){
			xCount[i] <- length(drop[drop >= x[i] & drop < x[i+1]])
		}
		for(i in 1:(length(x) - 1)){
			polygon(c(x[i], x[i+1], x[i+1], x[i]), c(-0.2, -0.2, -0.2+xCount[i]/max(xCount/0.25), -0.2+xCount[i]/max(xCount/0.25)), col = rgb(0, 0, 0, 0.1), border = FALSE)
		}
		axis(1)
		axis(2, line = -1)
		axis(4, at = c(-0.2, 0.05), labels = c('Min', 'Max'), line = -1)
		polygon(c(range, rev(range)), c(supHi, rev(supLo)), col = rgb(0, 0, 0, 0.25), border = FALSE)
		lines(range, supMed, col = rgb(0, 0, 0, 1), lwd = 2)
		mtext('Marginal Effect of Perceived Compromise on Incumbent Selection\nOver Ideological Distance from Alternatives', font = 2)
		mtext('Observed Values of Ideological Distance', side = 4)
		mtext('Change in Probability of Choosing Incumbent Party', side = 2, line = 1.5)
	dev.off()		

## calculate change in cabinet vote share ##

	data1 <- read.table('losses1.txt', sep = '|', header = TRUE, row.names = NULL)
	data2 <- read.table('losses2.txt', sep = '|', header = TRUE, row.names = NULL)
	data3 <- read.table('losses3.txt', sep = '|', header = TRUE, row.names = NULL)
	data4 <- read.table('losses4.txt', sep = '|', header = TRUE, row.names = NULL)
	data5 <- read.table('losses5.txt', sep = '|', header = TRUE, row.names = NULL)
	data6 <- read.table('losses6.txt', sep = '|', header = TRUE, row.names = NULL)
	data7 <- read.table('losses7.txt', sep = '|', header = TRUE, row.names = NULL)
	data8 <- read.table('losses8.txt', sep = '|', header = TRUE, row.names = NULL)
	data9 <- read.table('losses9.txt', sep = '|', header = TRUE, row.names = NULL)
	data10 <- read.table('losses10.txt', sep = '|', header = TRUE, row.names = NULL)

	data1$iteration <- paste(data1$k, '_1')
	data2$iteration <- paste(data2$k, '_2')
	data3$iteration <- paste(data3$k, '_3')
	data4$iteration <- paste(data4$k, '_4')
	data5$iteration <- paste(data5$k, '_5')
	data6$iteration <- paste(data6$k, '_6')
	data7$iteration <- paste(data7$k, '_7')
	data8$iteration <- paste(data8$k, '_8')
	data9$iteration <- paste(data9$k, '_9')
	data10$iteration <- paste(data10$k, '_10')
	
	data <- rbind(data1, data2, data3, data4, data6, data7, data8, data9, data10)
	votes <- c()
	for(i in 1:length(unique(data$iteration))){
		votes[i] <- mean(data$vote[data$iteration == unique(data$iteration)[i]])
	}

	summary(votes)
	summary(data)
	
pdf('totalEffects.pdf', width = 10)
	plot(1, 1, type = 'n',
		xlab = 'Predicted Incumbent Vote Change',
		ylab = 'Density of Predicted Effects',
		main = 'Predicted Incumbent Vote Change\nResulting from a One Standard Deviation Swing in Perceived Compromise',
		xlim = c(-0.1, 0.02),
		ylim = c(0, max(density(votes)$y)),
		axes = FALSE)
	axis(1)
	axis(2, at = c(0, max(density(votes)$y)), labels = c(0, 'Max'))
	polygon(density(votes), col = rgb(0, 0, 0, 0.25), border = FALSE)
	arrows(mean(votes), max(density(votes)$y)-3, mean(votes), 0, col = rgb(0, 0, 0, 1), length = 0.125, lwd = 2)
	text(mean(votes), max(density(votes)$y - 2), 'Mean', col = rgb(0, 0, 0, 1))

	arrows(quantile(votes, 0.96), (max(density(votes)$y)/3)-5, quantile(votes, 0.96), 0, col = rgb(0, 0, 0, 1), length = 0.125, lwd = 2)
	text(quantile(votes, 0.97), max(density(votes)$y)/5, '96th Percentile', col = rgb(0, 0, 0, 1))
dev.off()
